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We consider the general open system problem of a charged quantum oscillator confined in a 
harmonic trap, whose frequency can be arbitrarily modulated in time, that interacts with both an 
incoherent quantized (blackbody) radiation field and with an arbitrary coherent laser field. We 
assume that the oscillator is initially in thermodynamic equilibrium with its environment, a non- 
factorized initial density matrix of the system and the environment, and that ait = G the modulation 
of the frequency, the coupling to the incoherent and the coherent radiation are switched on. The 
subsequent dynamics, induced by the presence of the blackbody radiation and the laser field, is 
studied in the framework of the infiuence functional approach. This approach allows incorporating, 
in analytic closed formulae, the non-Markovian character of the oscillator-environment interaction at 
any temperature as well the non-Markovian character of the blackbody radiation and its zero-point 
fiuctuations. Expressions for the time evolution of the covariance matrix elements of the quantum 
fiuctuations and the reduced density-operator are obtained. 
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I. INTRODUCTION 

Since the seminal work of Magalinskii [1], Feynman 
and Vernon [2], UUersma [3], and Caldeira-Leggett [4], 
the theory of open quantum systems has been souree of 
great interest and object of continuous development, re- 
finement and application (cf. Refs. [5-9] and references 
therein). This theory has, for example, provided a solid 
conceptual framework to explain fundamental phenom- 
ena such as the quantum-classical transition [10], viola- 
tions of the second law of the thermodynamics [11], the 
survival of quantum features at high temperature [12, 13], 
and has found applications in several fields in physics and 
chemistry [5-9]. 

Despite the beauty and power of this theory, a study 
of the dynamics of a particular system can be very cum- 
bersome due to the sheer complexity of correctly incor- 
porating the various time, energy and coupling scales. In 
order to circumvent this problem various approximations, 
such as the weak-coupling-to-the-bath, Markovian, high 
temperature, or the initial factorizing condition is usually 
invoked [5-9, 14]. However, the development of modern 
experimental techniques for preparing and manipulating 
physical and chemical systems has reached the regime 
where such approximations are questionable. This fact 
has encouraged the development of techniques for more 
consistently treating and analyzing open quantum sys- 
tems (cf. Ref. [15]). 

Despite the refinement in technique, no approach is 
completely approximation- free [5-9, 15], and therefore 
some physical features of the dissipative dynamics arc of- 
ten absent in these descriptions [5-7, 15]. The breakdown 
of these common approximations is expected to occur in 
the low temperature regime [16] and/or in the presence 
of initial correlations between the system and the envi- 
ronment [14] and, in particular, in the case of driven 



non-equilibrium quantum systems [12, 17-19]. The ubiq- 
uitous presence of this situation (an open quantum sys- 
tems under the presence of time dependent fields in, e.g., 
coherent control scenarios of chemical systems [20] and 
assorted physical systems [19, 21, 22]), motivates a for- 
mal and detailed treatment of such dynamics. 

In this paper, we derive formal exact results for the 
non-Markovian dynamics of a prototypical system, in- 
cluding the presence of initial correlations between the 
system and the bath, and the possibility of arbitrary 
rapidly oscillating driving forces, and that can be ap- 
plied in the low temperature and strong coupling-to-the- 
bath regimes. The results derived here can be used to 
study, e.g., (i) the incoherent excitation of open quantum 
systems [21], a subject of great interest in the chemical 
physics community [23, 24]; (ii) environment-assisted one- 
photon phase control [22, 25]; (iii) optimal-control-based 
cooling of quantum nano-resonators by means of para- 
metric driving [19] and (iv) the establishment of a quan- 
tum limit on non-Markovian time scales in the spirit of 
Ref. [12]. 

The particular system studied here is a charged quan- 
tum oscillator confined in a harmonic trap that is initially 
in thermodynamic equilibrium with its environment (non- 
factorized initial conditions [14]). For i > 0, we start 
varying the frequency of the harmonic trap and couple 
the oscillator, via the dipole, with an incoherent quan- 
tized blackbody radiation field and with an arbitrary co- 
herent laser field. By means of the influence functional 
approach [2, 4, 14, 26], we derive analytic closed expres- 
sions for the non-Markovian time evolution of the covari- 
ance matrix elements of the quantum fluctuations valid 
at any temperature, any system-environment coupling 
strength and incorporating the zero-point fluctuations of 
the radiation. 
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II. DESCRIPTION OF THE SYSTEM 
HAMILTONIAN 

The Hamiltonian of a quantum oscillator confined in 
a harmonic trap and in contact with the environment 
comprises three parts: the Hamiltonian of the quantum 
oscillator itself, Hg', the Hamiltonian of the environment, 
Htb, here described according to the UUersma-Caldcira- 
Lcggett model [3, 4], as a collection of harmonic modes; 
and an interaction term _ffs_TB between the two systems. 
For typical bilinear coupling, the three contributions arc 
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with p and q the canonically conjugate momentum and 
position of the oscillator (an analogous notation describes 
the bath modes), m the mass of the quantum oscillator 
and a;^(i) = ujQ+ujp{t) the parametrically modulated fre- 
quency. This frequency comprises two components: wq, 
a constant frequency and an arbitrary time-dependent- 
frequency ujp{t). The magnitude of the interaction be- 
tween the system and the bath is determined by the cou- 
pling constants Cj. 

In the presence of the blackbody radiation the Hamil- 
tonian Hq = Hs + Hs-TB + needs to be augmented 
to include the interaction with the field as well as the 
Hamiltonian of the field modes, 

^BB = (^"LsOk.s + ■ (4) 

Assuming that the charged oscillator interacts weakly 
with each mode in the field, we can adopt a dipole-dipole 
type interaction, giving the overall Hamiltonian 

H^^(p-lA)\\m.itfe-qE,it) 

+ Hs^TB + Htb + Hbb, 

where e/c is the coupling constant to the radiation, ak,s 
and ^ are the annihilation and creation operators of 
the field mode of momentum k and polarization s. The 
vector potential is given by 

^ = E (^) ' /fc^k.., (flk,. + 4,.) . (6) 

where e is the polarization unit vector, V is the volume of 
the auxiliary cavity containing the field modes and is 
the electron form-factor (Fourier transform of the charge 



distribution) that incorporates the electron structure [27, 
28] . We have assumed, with no loss of generality, that the 
form factor and polarization vector are real. Note that 
by virtue of Eq. (4), Eq. (5) already contains zero-point 
or vacuum fluctuations. 

Equation (5) can be generalized to include an addi- 
tional term —qEi^(t), which allows for the possible ma- 
nipulation of the charged oscillator, via dipole coupling, 
by means of the electric field E'l (t) of a pulsed or contin- 
uous laser field. 

Since Eq. (5) includes the diamagnetic term A^, it 
is not suitable for a path integration calculation, which 
is why it is usually omitted [29, 30]. However, the 
contribution of this term is relevant for the derivation 
of the partition function of the oscillator in the pres- 
ence of the blackbody radiation (cf. the discussion in 
[31; 32]). In our case, this term can be introduced by 
means of the Power-Zienau's transformation [33] (see 
also [34] for the original version and [35] and references 
therein for a short historical review on this transforma- 
tion), T = exp |-^f <? • ^1 , which transforms p ^ p + |A, 

g -5- g, pk,s Pk,s + f mkWk9, gk,s "> gk,s, where we have 
defined [27, 28] 

a-k^s = (TOki^k^k.s + iPk,s)/v/2TOk^Wk, (7) 

with TOk = 47re^/^/(aJkV^)- The corresponding total 
Hamiltonian reads 

+ E^ + ^('^.-;^^j (8) 

+ E 0~~ ^P^'" + "^kWk?)^ + ^"Ik^k, Jk,s- 

'^k 

k,.s 

Here the oscillator is seen to be coupled to the momentum 
coordinate Pk,,s - 

Recalling that the electric field of the blackbody radi- 
ation is given by 

^=-a^ = iE( — ) Aek,.(ak,.-ay , (9) 

k.s 

and examining Eq. (7), we see that in Eq. (8), the posi- 
tion of the oscillator is coupled to the electric field of the 
blackbody radiation. From an open-quantum-systems 
perspective, this implies that any statistical behavior 
induced by the blackbody radiation (as seen elsewhere 
[21. 27; 28]) is dictated by the stochastic fluctuations of 
the electric field. 
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III. INITIAL DENSITY MATRIX 



After charactering the Haniihonian of the total system, 
the next step in the description of the dynamics is the 
determination of the initial state, taken to be the equi- 
librated state of S-I-TB. By denoting the coordinates of 
S-I-TB as Q = {q, Q}, the matrix elements of the initial 
density operator of the system S plus the environment 
TB can be calculated as [14] 

Q" 

(S" jp^l Q') - Z^^^ J 2?Qexp (-;^^^[Q]) , (10) 
C' 

where the integral is over all paths Q(t), < t < H/Stb 
with Q(0) = Q' and Q(;i/3tb) = Q" ■ The bar stands 
for the trajectories Q(t) in imaginary time t — > —ih/3, 
with /3tb = 1/(^81x8)- In the imaginary-time path 
integral described by Eq. (10), 5'^[Q] denotes the Eu- 

_ hf! _ ^ 

chdean action of the system S^[Q] = / drC^iQ, Q) = 



S'f [q] + 'S'f_xB Q] + "^TB [Q] ' obtained by introducing a 
global minus sign in the potential energy [14, 36]. Since 
Pl3 denotes the density operator at t = 0, we assume that 
a;p(0) = in Eq. (1). The matrix elements of initial 
density operator of interest are 

{qLQ" \m\ qLsQ') = {q'L \h.J<lLs) {Q" \pp\ Q') , 

(11) 

where p^g^ denotes the equilibrium density operator of 
the cavity only, at temperature Tbb and defined as in 
Eq. (10) using the Euclidean action ^^^[gk.s] of the bare 
cavity, pp is the thermal density operator of system and 
bath introduced in Eq. (10), and the blackbody-radiation- 
mode-coordinates gk,s are defined in Eq. (7). When a 
system S is in contact with two thermal baths, what is 
usual in the literature (cf. [21, 22]) is that the initial state 
of the total system is assumed to be factorized, in our case 
this is equivalent to take p(0) = PPbb ® PPtb '■^ PS, being 
ppT-B and ps the density operator of the thermal bath 
TB and the system S, respectively. Here we deal with a 
more complex situation because we take into account the 
initial correlations between the system S and the thermal 
bath TB. 

In general, one would also like to study the system S 
prepared in a state other than the equilibrium state pp, 
e.g., in a coherent state or in a squeezed state. Accord- 
ing to Refs. [14, 37], one could prepare a different ini- 
tial state by allowing the operators An , acting only 
in the system Hilbert space of S to generate an ini- 
tial non-equilibrium density operator of system and bath 
Pa = J2ni^n iTB)p/3(iTB <E) A'^) , where Itb denotes 
the unit operator in the Hilbert space of the bath alone. 
In the position representation, the matrix elements of pp 



are given by 

\p',\ Q'_) = J dQ"dQ'X{q'^,q",q'_,q') 

xd{Q'+-Q")5{Q'_-Q'){Q"\pp\Q'), 

(12) 

where the propagating function X{q'^,q" ,q'_,q') = 

J2n{l'+\^n\q"){q'\A',^\q'-) characterizes the action of 
these operators. The delta functions indicate that the 
imaginary-time paths for the bath degrees of freedom are 
continuously connected to the real-time paths describing 
the time evolution of the initial state [14, 37]. The ther- 
mal initial state in Eq. (10) can be recovered by setting 
the auxiliary operators An and A' to Is, which yields 
A(g+, g", qL,q') = S{q'_^_ - q")S{q'_ - q'). As distinct from 
Q' and Q" , which are the endpoints of a single imaginary- 
time-trajectory, Q'_^ and Q'_ denote the initial condition 
for two real-time-trajectories and "— " (see below). 
For the case of unitary evolution these can be identified 
as the forward and backward trajectories associate to the 
unitary time-evolution-operator and its adjoint, respec- 
tively. Additional details about different initial prepara- 
tions can be found in Refs. [14, 37]. 



IV. DYNAMICS OF THE SYSTEM 



In this section we solve for the time evolution of the 
initial density matrix [Eq. (11)] under the action of the 
Hamiltonian [Eq. (8)] , using the Feynman and Vernon in- 
fluence functional approach [2, 26]. For this problem we 
require a mixture of the influence functional for factor- 
izing initial conditions [2, 4, 17, 26] in order to consider 
the effect of the radiation, and for non-factorizing initial 
conditions [14, 37] in order to correctly characterize the 
equilibrium state between the oscillator and the thermal 
bath at t = 0. 

The standard path integral calculations are performed 
for position-position couplings [2, 4, 5, 14, 17, 26, 37] 
[cf. the term qqj in Eq. (8)], and for momentum- 
momentum coupling [4]. The oscillator-radiation cou- 
pling is of position-momentum type, qpk.s [cf. Eq. (8)]. 
By means of a set of unitary transformations [28], one 
could invert the role of the momentum and position op- 
erators of the field [see Eq. 7], with the caveat that this 
generates an additional term proportional to the initial 
condition of the oscillator (see below). However, since 
the system described in Eq. (8) is still linear, the path in- 
tegral calculation can also be carried out analytically for 
the position-momentum coupling. This is the approach 
followed below. 
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A. Derivation of tlie propagating function and the 
influence functional 



The time evolution of the system S is described by the 
reduced density operator ps(t) — trTB,BBp(^)- FoUowing 
Ref. [14], we obtain that the matrix elements of /5s(t) are 
given by 



(<Z+|ps(i)k-) = / dqVdg'_dg'WA(g;,g",g^,g') 
X J(<z:^,g^,t;g;,g^,0;q-",g'), 



(13) 



where J((j" , , i; g!., 0; g", g') is the propagating func- 
tion of the system density matrix which can be expressed 
in terms of the functional phase S [(7+ , <Z- , by means of 
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the three-fold path integral expression 
J{q'l,q'L,Uq'^,q'_,0;q'\q') = 
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where Z normalizes J {q'_l , q'_L , t; q'_^,q'_,0; g", q') to 5{q'[_ — 
q'+)5{q'L-q'_){q"\pp,sW) at f = 0, being pp^s =trTB(/0;3). 
The real time path-integrals over q^ and q- are carried 
out subject to the endpoints (?+(0) = g^, q+{t) = q+, 
g_ (0) = q'_ and g_ {t) = g" , while the imaginary time 
path integral are over g(0) = q' and ^(/i/^tb) = q" ■ Re- 
call that the imaginary time path-integral allows for the 
calculation of the equilibrated density operator of S+TB 
and the influence of their initial correlations in the sub- 
sequence time evolution. 

After tracing over the degree of freedom of TB and BB, 
and after defining q^ ~ [r + x)/2 and g_ ~ r — x, we 
have that the functional phase S[a;, r, q] is given by 



dakTBir ~ a)q{T)q{a) 



+ dr dsK^jiis - ir)g(r)a;(s) 



ds \^inx{s)'r{s) — mujQr{s)x{s) — mujp(t)'^r{s)x{s) + Ei^(t)x{s) — r'riTB{s)x{s)} 

t 

du [?7tb(s) + r/BB(s)] x{s)r{u) ^ J ^'^ [^''^tb(* - u) + ^bb(s - u)] x{s)x{u) 



(15) 




Note that r' = q"^ + ^q'_ , r" = q'i^ + \q'L and analogously 
for x" and x' . The various kernels entering into Eq. (15) 
are defined in the next section. 

The first term in the first line of Eq. (15) accounts for 
the equilibrium density operator of S in the presence of 
the thermal bath TB while the second term containing 
q{T)x{s) is responsible of the the effect of initial corre- 
lations between the environment and the system on the 
subsequent time evolution. The first three terms in the 
second line are responsible of the evolution under the 
parametric harmonic potential, while the fourth term is 
responsible for the evolution induced by the laser field E'l ■ 
As previously noted, we take a;p(0) = and -E'l(O) = 0. 
The last term in the second line of Eq. (15) arises from the 
incoherent excitation induced by the position-position 
coupling to the thermal bath TB. Since the coupling to 
the blackbody radiation field BB is of a different nature, 
position-momentum coupling, this transient term propor- 
tional to the initial position r' is not present. However, if 
one changes the role of position and momentum, as dis- 



cussed above, this transient term enters implicitly. The 
terms in the third line constitute the exponent of the in- 
fluence functional of the Fcynman- Vernon theory under 
the action of the thermal bath TB and blackbody radia- 
tion BB. 

The additional time integration in the last line of 
Eq. (15) accounts for the non-local time (non-Markovian) 
evolution of the density operator. Although the tempo- 
ral non-locality is determined by the various kernels in 
a cumbersome way, we can identify two kinds of non- 
Markovian contributions: one from the dissipative part 
and determined by the non-local character of 77TB and 
77BB, and a second determined by the thermal fluctua- 
tions described (see below) by the kernels and /^bb ■ 
The presence of the latter is not determined by the pres- 
ence of the former, i.e., in the limit of local dissipative 
Ohmic kernels, 77TB ~ ^(s), the non-local character of 
the thermal fluctuations is still present; it only vanishes 
in the high temperature regime [5, 38]. 
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B. Kernels in the functional action 

The quantities introduced in the effective action 
E[x, r, ga:] are defined in terms of the bath spectral density 
Jtb and the blackbody-radiation spectral density Jbb- 
These spectral densities are determined [4, 14, 27, 28] 
from the parameters of the bath modes and the coupling 
constants by means of 



OO 



Jbb(w) = ^'^m]^u!lS{u! - cjk) 



(16) 



(17) 



k,s 



Assuming that the thermal bath is dense in the frequency 
of the modes [4] , it is customary to describe the spectral 
density in Eq. (16) by, e.g., assuming the Ohmic model 



Jtb(w) = m7TBWf^TB/(^TB + ^^)) 



(18) 



where 7tb is the energy decay rate and ^Itb is a cutoff 
parameter related to the inverse of the bath memory time. 
In contrast with the thermal bath case, no assumption on 
the functional form of the spectral density of the black- 
body is needed in the continuous limit [27, 28, 30, 39]. 

The transvcrsality condition implies that only two of 
the three components of k contribute to the coupling 
[28], giving a global factor of two-thirds for the spec- 
tral density in Eq. (17). In the continuous limit, ~^ 
J dk, the spectral density for the blackbody radia- 



tion is 



>/bb(w) = Mtbb^^^^bb/ (^bb 



(19) 



where M ~ m + Mtbb^bb is the rcnormalizcd mass, 
tbb = 2e^/3Mc'^ and J^bb is a frequency cutoff. This 
path-integral-based expression coincides completely with 
the seminal results in Refs. [27, 28, 39] using the quan- 
tum Langevin formalism. It also coincides with the re- 
sult derived in Ref. [30] using the path integral approach. 
However, we need to note that in [30], the system, an 
electron, is interacting with its own radiation; here, by 
difference, we consider the system as being irradiated 
by an external blackbody radiation such as sunlight or 
moonlight for t > 0. This is precisely what allows us to 
separate the initial density operator of the system and 
the radiation. This natural emerging functional form of 
Jbb(i^) reveals, from a statistical viewpoint, the intrinsic 
non-Markovian character of the radiation [27, 28, 30, 39]. 
This fact implies that the two point correlation function, 
{E{t")E{t"))BB: of the electric field in Eq. (9) is not delta 
correlated. From an optics point of view [40] this means 
that the blackbody radiation is coherent, although the 
coherence time is very short, ~8 fs at Tbb = 5900 K (cf. 
Ref. [21] and Chap. 13 in Ref. [41]). 

Once we have condensed the relevant information for 
the thermal bath and the radiation field in the spectral 



densities Eqs. (18) and (19), we are in the position to 
define the various functions entering into the functional 
actions Yi[x,r^q]. 

The kernels _fCTB(s — it) and A:tb(''') arc given by [14] 

i^TB (s - it) = IC^B (s - ir) + iKi^B - ir) , (20) 



fcTB(T) 



TB 



Cn(0)exp(i:.„T) (21) 



with the Matsubara frequencies Un = 27rn//i/3TB- The 
real and imaginary parts of K^b{s — It) are given by 

OO 

j^rc . /■ 7 ^ cosh[a;(in/3TB) - t] 

^tb(s - 1^) = / — Jtb(w) . ■ ,\ ■ ^ , — cos(ws) 

J TT smh(ia;ft^TB) 

OO 

m ^—^ 

5„(s)exp(ii/„r), (22) 



TB 



do; sinh[aj(ift/3TB) — t] . 



KtUs - It) = - / — Jtb(^) ? /.V^r ^ ' sin('^s) 
TT smh(ia;ftpTB) 




m 



TB 



i/„(s)exp(ii/„T). 



(23) 



n— — OO 



while Cn{s) = ± r duJr^M ^^cos(tjs) or in 



terms of the damping kernel 7tb(=^) [14], Cn{s) — 

OO 

i|i/„| / du7TB(w)[exp(-|i^s(s + u)\) + cxp(-|i^s(s - u)\)]. 



The functions (?n(s) and /n(s) can be expressed in terms 
of the damping kernels 



2 fdu 
7tb,bb(s) = — / — 
m J TT 



2 f duj Jtb,bb(w) 



cos(ws), (24) 



and Cn(s) as 5„(s) = 7tb(s) - Cn(s) and /„ = 

The spectral density in Eq. (18) generates the dis- 
sipativc kernel 7(5) = 7TBf^TB cxp(— riTB|s|)- In the 
limit when the cutoff frequency ^Itb tends to infin- 
ity, 7(5) — ^ 27tb^(s), which corresponds to Markovian 
Ohmic dissipation. The spectral density in Eq. (19) gen- 
erates 7bb(s) = tbbS^Ib [2'5(s) - f^BBexp(-riBB|s|)] . 
Note that there is a fundamental limitation to the use 
of Eq. (19). That is, in the limit IIbb ^ 00, we get 
the surprising result that 7bb(s) = 0, i.e. no relaxation 
[42, 43]. This corresponds to the point-electron limit 
[fk = f^BB/(^BB + ^fe) = 1 in Eq. (6)] and is unphys- 
ical because even for the electron, 51bb remains finite, 
although large. According to Refs. [42, 43], there is a 
natural upper value given by ^Ibb ~ tbb^^, which corre- 
sponds to two-thirds of the time for a photon to traverse 
the classical electron radius (r^j = 2.818 x 10~"'^^m). Be- 
yond this natural limit, causality is violated [42] and the 
bare mass m takes negative values [42]. 
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Finally, the kernel Kbb{s) is given by 

oo 

-ft^BB(s) = / — JBB(w)cotll 



cos(ws). (25) 



The real parts of the kernel K^Bis — It) and K^b^s) are 
responsible for the decoherence due to thermal fluctua- 
tions, while the kernels 



lution to the equation of motion 



mf + TO— / du7(s — it)r(w) + TOaj(t) r 
ds 



??tb(s) m7TB(s) 
'?bb(s) = »7i,7bb(s) 

result in the relaxation process. 



(26) 
(27) 



mi — TO— / du7(s — u)2;(m) + 7Ticj(t)^a; = 0, 
ds ' 



(31) 



(32) 



C. Explicit calculation of the propagating function 



The explicit calculation of the propagating function 
demands evaluating the path integral in Eq. (14). Since, 
the system is linear, the path integral can be performed 
by evaluating the action in Eq. (15) along its stationary 
trajectories and condensing the effect of the fluctuations 
in a global time dependent factor [5, 14]. The extrema 
of the action are given by 



where we have defined £'l(s) = E\^{s) + ?''77bb(s) + 
ra\xC\(s)—\xC-2,(s')\ withf = (g-|-q')/2 and a; = q'. Ad- 
ditionally, we have defined r\[s) = mj{s) and R{s,u) = 
Rtb{s, u) + K{s — u)/m with 

i?TB(s,w) = -AtbCi(s)Ci(w) 



+ 



fcfl XI ""[5«(s)ff«(w) - /„(s)/„(ti)]. 



Ci{s) - 
C2{s) = 



1 



S/StbAtb 
1 



n— — oo 



h(3- 



TB 



(33) 

(34) 
(35) 



n— — OQ 



mq—mujQq- 



where Atb = 7^ J2n=~oo and u„ = (w^ + f2+C„) ^ 

dcrfcTB(T-cr)g(cr) = -i / dsKi^{s~iT)x{s) , gj^^^,^^ ^ harmonic potential, the functional action 

r, g^^] can be evaluated using only the real part of the 

(28) 

trajectories r{s) and x{s) [14, 37, 44, 45], we need to solve 
only for the real part of Eq. (31) and (32). Due to the 
linear character of (31), the solution to the homogeneous 
part can be written as 



mr + mcuCt) r — Ei^(t)r -\ / durjis — u)r(u) 

as ' 



f'f;BB(s)+i / duK{s - u)x{u) + 



(29) 



r{s) 



Mt) 



+ r'[Ms) 



Ms) , (36) 



where (/)i(s) is the fundamental solution for r(0) = and 
r(0) ~ 1, while (j)2{s) is the fundamental solution for 
r(0) = 1 and r(0) = 0. Thus, for Eqs. (31) and (32) we 
have 



mx + muj{t)'^x - 



d_ 

dl 



duri{s — u)x{u) — 0, (30) 



where we have defined r]{s) ~ ?yTB(s) + '7bb(s) and 
K{s) = K^q{s) + A'bb(s). The term r'ryBB(s) appears 
here as a consequence of the sudden turn-on of the black- 
body radiation. Since, we assume that the parametric 
driving, as well as the laser field and the blackbody radi- 
ation act after t — 0, the equilibrium state of our system 
coincides with the one derived in Ref. [14] , so we need to 



+ — / duMs - u)F'{u) 
m Jq 

1 Ms) 



^it) 



lit) 



Ms) 



x{s) = X 



m 01 (i) Jo 
,fMs) 



du(f>i{t - u)F'{u), 



Mt) 



+ x'[Ms)-^,Ms) 

Mt) 



(37) 



(38) 



focus only on the evaluation of the real part of the action. , ^ \ • r 1 ^ 1 1 i- r fn\ n i 

where 1^31(5) is the fundamental solution tor xiv) = (J and 

The real part of the action is stationary along the so- i(0) = 1, while 1^32(5) is the fundamental solution for 
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x(0) = 1 and a;(0) = 0. For the particular case of no 
frequency modulation, ^i(s), <p2{s), ipi{s) and 1^2(5) can 
be derived from standard Laplace techniques [14]. For 
Markovian dissipation, wtb 00, and harmonic mod- 
ulation of the frequency, 4>i{s), 4>2{s), 'Pi{s) and 'fi2{s) 
are related to the Mathieu functions [12, 17]. For more 
general cases, these functions must be calculated numeri- 
cally. However, the functional form of Eqs. (37) and (38) 
is very convenient of the subsequent analytical calcula- 



tions. For further convenience we define 
vi[t,s) = ip2[s) — (pi(s), V2[t,s) = — --, (39) 



ui(i,s) = 02(s)-^0i(s), ^2(i,s) = Mfl. (40) 



01 w 



v{t) 



The influence functional in Eq. (14) can now be rewrit- 
ten as 



J{r",x",t;r',x',Q]f,x) 



(41) 



where Nit) is a normalization factor given by N{t) = 

2T:h^\u2{t,Q)\{2'Kh^KY''^ . After evaluating Eq. (15) 
along r''°(s) and x{s), we get 



v^r " // / / - -1 

L[r , a; , r , a; ,r, xj = im 



-2 ^ i::TB ^2 j _^ rn{x"r"u2{t, t) + x'r'ui(t, 0)) - m {x'r"u2{t, 0) - x"r'ui{t, t)) 



2A 



TB 



-Ei^{s) + r'7BB(s) + fCi{s) - \xC2(,s) 



+ m J ds [x'vi (t, s) + x"v2 (t, s)] 


+ —mx'^ ds duR{s,u)vi(t, s)vi{t,u) + —mx'x" ds duR{s,u)vi{t, s)v2{t,u) 
2 Jo Jo 2 Jo Jo 

+ —mx"x' ds duR{s,u)v2{t, s)vi{t,u) + —mx"^ ds duR{s,u)v2{t, s)v2{t,u), 
Jo Jo 2 Jo Jo 



(42) 



■ Si^-oo '""('^0 + Cn)- These results no modulation of the frequency and for factorized initial 



where Qtb = 

provide the general expression for the influence func- conditions, Eq. (41) reduces to [4, 48] (see also [45] for a 



tional. 



description in terms of the Wigner function) 



D. Limiting cases 

The general result in Eq. (41) includes, and agrees 
with, several limiting cases. These include: 

— In absence of blackbody radiation and with no para- 
metric modulation of the frequency, Eq. (41) reduces to 
the result in [14]. 

— In absence of the thermal bath and for no laser field 
and no parametric modulation of the frequency, Eq. (41) 
is the formal path integral equivalent of [27, 28] (based 
on the quantum Langcvin equation formalism) . 

— In absence of the thermal bath and with no paramet- 
ric modulation of the frequency, Eq. (41) is the formal 
path integral equivalent of [27, 28] to the result in [46, 47]. 

— In absence of blackbody radiation, for no laser field, 
harmonic modulation of the frequency and factorized ini- 
tial conditions, Eq. (41) reduces to [17] (see also [12]). 

— In absence of blackbody radiation, for no laser field, 



E. Evolution of the density matrix 

For an initial thermal state, i.e. X{q'_^,q" ,q'_,q') = 
5{q'j^—q")5{q'_—q') in Eq. (12) and correspondingly f = r' 
and X = x' , the influence functional in Eq. (41) with the 
function phase given in Eq. (42) can be written in the 
very compact form 



J{r",x",t;,r',x',0;r',x') 



Nit) 



X exp |-ix'^A(t)x + El{s)x' + E^{s)x"^ 



(43) 



where A{t)K = xjA/(t)x/ -I- x^Ai(t)xi -I- 2Bxi where 
X = ix",r",x',r'), X, = (x',r'), x/ = {x",r"), B = 
(A13X" -|- A23r", Ai4a;") with the time dependent matrix 
A{t) given by 
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( i?22(i) -iu2(i,i) -C^{t) + Ri2(t) -mi(i,t) -iCf (i) \ 

-\U2{t,t) iw2(i,0) 

-C^ {t) + Ri2{t) m2(i,0) Otb -2C2+W+i?ii(t) -mi(t, 0) - i(7+(t) 

V -mi(t,t) -iCf (t) ~mi{t,Q)-\C+{t) 1/Atb / 



(44) 



being 

t f 
C+(t)= y AsC,{s)v,{t,s), CJit)^ 





dsCj(s)t;2(i, s), 
(45) 



t t 

-B^(i) = j dsEi^{s)vi{t,s), E~{t) = J dsEL{s)v2{t, s), 



(46) 



i?y (t) 



ds / duR{s,u)vi{t, s)vj{t, s). 



(47) 



Cf*^ (t) is obtained by replacing Ci (i) — > Ci (i) + 7bb (i) in 
the definition of C^it) in Eq. (45). 



r 



The first moments are determined by 



1 



(pit)) = / ds 

1 



mu2{t,0) Jq 
t 



dsvi{t, s)£'l(s) 



(52) 







"2(i,0) Jo 



Mt,s)^^^^+V2it,s) 
U2[t,0) 



dsvi{t,s)Ei^{s), 



(53) 



from Eqs. (52) and (53), one can verify that {p{t)) 
md{q{t))/dt. 



The second moments arc given by 



ft2 



Mn(t)--Mi2(i)' 



(54) 



V. TIME EVOLUTION OF THE COVARIANCE 
MATRIX ELEMENTS 



Due to the linear nature of the system under considera- 
tion, the dynamics as well as the statistical properties can 
be characterized in terms of the first and second moments 
[48], defined as {f{q{t))) ^ J dr" fir")ps{r" , x" = 0,t), 
or in terms of the propagating function in Eq. (43) by, 

if (lit))) ^ J dr"dr'dx'f{r")J{r'\x" = 0, r', x', r', x')- 

(48) 

So, we can now make use of the explicit form of the prop- 
agating function in Eq. (43). Note that the propagating 
function in Eq. (43) also contains the initial thermal state 
Pis^s- These moments define the variances or dispersion 
relations 



(49) 



^,pit) = T^mvit) +pit)q{t)) - m){p{t)), (50) 



(51) 



which will be used in order to express explicitly the time 
evolution of the density matrix in Eq. (58). 



-{pq + qp) =m{q^{t))u2{t,t) 



U2{t,Q) 
. ft^ATB 



■mu2{t, 0) 



[R,2{t)-C^{t)] 

Cf(t)+ui(t,d Mi2{t) 



(55) 



+ {q{t)) / dsv2{t,s)Ei^{s), 
Jo 

{p^{t)} = hmR22{t) + /imAxB (t) + ui{t,t) 
- m^iilit, t){q^{t)) + mu2{t, 0){pq + qp) 

+ / dsv2{t,s)Ei^{s) , 
Jo 

where 



(56) 



'TB 



2C+{t) + Rii -i ui(t, 0) + C+(t) 



(57) 



— 1 



?ii(t,0) + C+(t) 



l/A- 



TB 



In the absence of the parametric driving, the black- 
body radiation and the laser field, Eqs. (54)- (56) are time- 
independent and coincide with the expressions (6.62)- 
(6.64) inRcf. [14], i.e., (q^it)} = (?i/to)Atb, {pq + qp) = 
and {p'^{t)) = hmilTB- 
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VI. TIME EVOLUTION OF THE REDUCED 
DENSITY-OPERATOR 

In terms of the second moments, the evolution of the 
reduced density operator can be written as 



cxp 



1 



<Tpp{t) 



„"2 



[r"-{qm 



(58) 

In summary, to obtain the time evolution of the reduced 
density operator, we proceed as follows: 

(i) One first specifies the spectral density Jtb(w) to 
obtain the function describing the modulating force in 
Eq. (8). This permits us to obtain the fundamental so- 
lutions 01. 2(5) and (y9i^2(s) in Eq- (37) and Eq. (38), re- 
spectively. 

(ii) One then calculates the kernels Ktb{s), k{T) and 
Kbb{s) defined in Eqs. (20), (21) and (25). 

(iii) With the fundamental solutions obtained, and all 
the kernels calculated, we calculate the auxiliary func- 
tions wi,2(^, s) and ui^2it,s) given in Eqs. (39) and (39), 
and subsequently the functions C^{t), E^{t) and Rij(t) 
defined by Eqs. (45), (46) and (47), respectively. 

(iv) One then calculates the first and second moments 
given in Eqs. (52-56), and subsequently the dispersion 
relations in Eqs. (49-51) and system dynamics via Eq. 
(58). 

This brief prescription concludes our completely formal 
and approximation- free treatment. Note that we have 
successfully applied the method to a number of cases, 
some of which will be reported elsewhere [21, 22]. 



insight into the role of coherences detected in photosyn- 
thetic light-harvesting complexes (for a review in the sub- 
ject see [15]). In doing so, we need to translate the propa- 
gating function in Eq. (43) into the energy basis in order 
to identify the incoherent/coherent nature of the excita- 
tion [21]. 

Our results allow us to directly study the possibility of 
environmentally assisted one-photon phase control [25] 
provided by the fact that the initial equilibrium density 
matrix deviates from the canonical distribution. In this 
respect, we eliminate the incoherent radiation and the fre- 
quency modulation contributions and focus on how the 
phase information encoded in Ei^(t) can be used to ma- 
nipulate the populations of the oscillator [22]. 

Additionally, the analytic closed expression could be 
useful in understanding the delicate balance between dis- 
sipation and driving under non-Markovian evolution that 
has been pointed out in Rcf. [19] in the context of optimal 
control theory and cooling of nano-mechanical resonators. 
In particular, the optimal cooling protocol addressed in 
Ref. [19] by means of numerical techniques, can be ana- 
lyzed in great detail from the second moments derived in 
Eqs. (54)-(56) and the theory of variational calculus. 

In Rcf. [12], it was established that the usual quantum 
limit, huj/kBT > 1, needs to be reformulated for out-of- 
equilibrium systems. However, in that work a Marko- 
vian Ohmic spectral density, r^xB cxd in Eq. (18), was 
used, disregarding in this way the dynamics during non- 
Markovian time scales. The physical system considered 
in [12] consisted of two identical harmonic oscillators with 
time-depend coupling, in the normal mode description 
we get two independent parametric oscillators. So, under 
the same circumstances considered in [12], the results 
derived here allow us to explore the limit for the pres- 
ence of quantum features in non-Markovian-driven-open- 
quantum systems, which is of great importance in, e.g., 
quantum statistical mechanics or control theory. 

Studies of these applications are underway in our lab- 
oratories. 
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